Load required packages for entire analysis
dlnm.pkg <- c("tidyverse", "lubridate","xts","ggthemes","PerformanceAnalytics","reshape2","rugarch","timetk","parallel","timeSeries","tseries","data.table","ggplot2","dlnm","broom","caret","gridExtra","splines","splines2","pspline","cowplot","mgcv","spi","chron","gridGraphics","grid","pscl","MASS", "AER", "Hmisc", "MuMIn", "VGAM", "forecast", "seasonal", "plotly", "ggmap", "rgeos", "tmap", "maptools", "maps", "ggfortify", "htmltools","webshot","knitr","flexdashboard", "imager", "httr", "curl", "gmodels","jpeg","rgdal","raster","sp","here")
lapply(dlnm.pkg, library, character.only=TRUE)
Load Malawi map dataset
dlnmtmp <- tempfile()
download.file("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/malawi_map.zip", destfile=dlnmtmp)
unzip(dlnmtmp, exdir = ".")
malawi.map <- rgdal::readOGR(".","malawi_map")
Subset Blantyre map from Malawi map
blantyre1.map <- malawi.map@data$OBJECTID >289 & malawi.map@data$OBJECTID <297
blantyre2.map <- malawi.map@data$OBJECTID >308 & malawi.map@data$OBJECTID <311
blantyre3.map <- malawi.map@data$OBJECTID >342
Convert shape files into dataframe
blantyre.map <- rbind(fortify(malawi.map[blantyre1.map,]), fortify(malawi.map[blantyre2.map,]), fortify(malawi.map[blantyre3.map,]))
blantyre.map$id <- as.integer(blantyre.map$id)
Load Blantyre demographics’ and features’ dataset
blantyre.demog <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_demog.csv"))
map.features <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/blantyre_features.csv"))
blantyre.demog$id <- as.integer(blantyre.demog$id)
map.all <- merge(x=blantyre.map, y=blantyre.demog, by="id", x.all=TRUE)
rm(list = ls()[grep("^blantyre", ls())])
Load case dataset and subset by serovar
case <-read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/case.csv"))
case$case_date <- dmy(case$case_date)
case$case_count <- c(1)
case.typhi <- subset(case, organism=="typhi")
case.iNTS <- subset(case, organism=="iNTS")
Assign 0 when no case was diagnosed on that day
case.typhi <-aggregate(case.typhi$case_count, by=list(case.typhi$case_date), FUN=sum, na.rm=TRUE)
names(case.typhi) <- c("date", "case_count")
case.typhi <- merge(case.typhi, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.typhi[is.na(case.typhi)]<-0
case.iNTS <-aggregate(case.iNTS$case_count, by=list(case.iNTS$case_date), FUN=sum, na.rm=TRUE)
names(case.iNTS) <- c("date", "case_count")
case.iNTS <- merge(case.iNTS, data.table(date=seq.Date(min(case$case_date), max(case$case_date), by="day")), by="date", all=TRUE)
case.iNTS[is.na(case.iNTS)]<-0
Convert to XTS and aggregate monthly
case.typhi = as.xts(case.typhi[,-1,drop=FALSE], order.by=as.Date(case.typhi[,1]))
case.iNTS = as.xts(case.iNTS[,-1,drop=FALSE], order.by=as.Date(case.iNTS[,1]))
case.typhi <- apply.monthly(case.typhi, FUN=sum)
case.iNTS <- apply.monthly(case.iNTS, FUN=sum)
Load Blantyre climate dataset
climate <- read.csv(curl("https://raw.githubusercontent.com/deusthindwa/dlnm.typhoid.nts.climate.blantyre.malawi/master/data/climate.csv"))
Define average rain/temp between chileka & chichiri stations
climate$climate_date <- dmy(climate$date)
climate$rainfall <- (climate$chil_r+climate$chic_r)/2
climate$temperature <- ((climate$chil_mint+climate$chil_maxt)/2+(climate$chic_mint+climate$chic_maxt)/2)/2
climate <- subset(climate, select=c(climate_date, rainfall, temperature))
climate.rain <- subset(climate, select=c(climate_date, rainfall))
climate.temp <- subset(climate, select=c(climate_date, temperature))
Convert to XTS and aggregate monthly
climate.rain <- as.xts(climate.rain[,-1, drop=FALSE], order.by=as.Date(climate.rain[,1]))
climate.temp <- as.xts(climate.temp[,-1, drop=FALSE], order.by=as.Date(climate.temp[,1]))
climate.rain <- apply.monthly(climate.rain, FUN=mean)
climate.temp <- apply.monthly(climate.temp, FUN=mean)
Create datasets for time series decomposition
case.typhi <- tk_tbl(case.typhi, preserve_index=TRUE, rename_index="date")
case.iNTS <- tk_tbl(case.iNTS, preserve_index=TRUE, rename_index="date")
climate.rain <- tk_tbl(climate.rain, preserve_index=TRUE, rename_index="date")
climate.temp <- tk_tbl(climate.temp, preserve_index=TRUE, rename_index="date")
Seasonal-adjusted NTS cases
case.iNTS.ts <- ts(na.omit(case.iNTS$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.iNTS.ts)))), preserve_index=FALSE)
case.iNTS$case_count_sea <- trend_n$value
case.iNTS$case_count_sea[case.iNTS$case_count_sea < 0] <- 0
setnames(case.iNTS, old="case_count", new="case_count_obs")
trend_n <-tk_tbl(exp(mstl(log(case.iNTS.ts))))
case.iNTS$case_count_tre <- trend_n$Trend
Seasonal-adjusted typhoid cases
case.typhi.ts <- ts(na.omit(case.typhi$case_count), frequency=12)
trend_n <- tk_tbl(exp(seasadj(mstl(log(case.typhi.ts+1)))), preserve_index=FALSE)
case.typhi$case_count_sea <-trend_n$value
case.typhi$case_count_sea[case.typhi$case_count_sea < 0] <- 0
setnames(case.typhi, old="case_count", new="case_count_obs")
trend_n <- tk_tbl(exp(mstl(log(case.typhi.ts+1))))
case.typhi$case_count_tre <- trend_n$Trend
Seasonal-adjusted rainfall
climate.rain.ts <- ts(na.omit(climate.rain$rainfall), frequency=12)
trend_n <- tk_tbl(seasadj(mstl(climate.rain.ts)), preserve_index=FALSE)
climate.rain$rainfall_sea <- trend_n$value
climate.rain$rainfall_sea[climate.rain$rainfall_sea < 0] <- 0
setnames(climate.rain, old="rainfall", new="rainfall_obs")
trend_n <-tk_tbl(mstl(climate.rain.ts))
climate.rain$rain_tre <- trend_n$Trend
Seasonal-adjusted temperature
climate.temp.ts <- ts(na.omit(climate.temp$temperature), frequency=12)
trend_n <-tk_tbl(seasadj(mstl(climate.temp.ts)), preserve_index=FALSE)
climate.temp$temperature_sea <- trend_n$value
climate.temp$temperature_sea[climate.temp$temperature_sea < 0] <- 0
setnames(climate.temp, old="temperature", new="temperature_obs")
trend_n <- tk_tbl(mstl(climate.temp.ts))
climate.temp$temp_tre <- trend_n$Trend
rm(list=ls()[grep("^trend_n", ls())])
Plots of seasonal-adjusted incidences (Supplementary Figure 2)
Th <- theme(plot.title=element_text(hjust=0)) +
theme(axis.title.x=element_text(size=10)) +
theme(axis.title.y=element_text(size=10)) +
theme(axis.text.x=element_text(face="bold", size=10), axis.text.y=element_text(face="bold", size=10)) +
theme(legend.key.height=unit(1, "line")) +
theme(legend.key.width=unit(1, "line"))
S1 <- ggplot(as.data.frame(case.iNTS)) +
geom_line(aes(date, case_count_obs, color="Original"), size=0.8) +
geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="orange2")) +
labs(title="A", x="Year", y="NTS cases") + Th +
theme(legend.justification=c(0.5,0), legend.position=c(0.7,0.7), legend.text=element_text(size=10),
legend.title=element_text(size=0))
S2<-ggplot(as.data.frame(case.typhi)) +
geom_line(aes(date, case_count_obs, color="Original"), size=0.8) +
geom_line(aes(date, case_count_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="red2")) +
labs(title="B", x ="Year", y="Typhoid cases") + Th +
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10),
legend.title = element_text(size=0))
S3<-ggplot(as.data.frame(climate.rain)) +
geom_line(aes(date, rainfall_obs, color="Original"), size=0.8) +
geom_line(aes(date, rainfall_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="blue2")) +
labs(title="C", x ="Year", y = "Rainfall (mm)") + Th +
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10),
legend.title = element_text(size=0))
S4<-ggplot(as.data.frame(climate.temp)) +
geom_line(aes(date, temperature_obs, color="Original"), size=0.8) +
geom_line(aes(date, temperature_sea, color="Seasonal-adjusted"), size=0.8) +
scale_color_manual(values = c("Original"="black","Seasonal-adjusted"="green2")) +
labs(title="D", x ="Year", y = "Temperature (°C)") + Th +
theme(legend.justification=c(0.5,0), legend.position = c(0.3,0.7), legend.text = element_text(size = 10),
legend.title = element_text(size=0))
grid.arrange(grobs=list(S1, S2), ncol=2, nrow=1)

Calculate NTS/typhoid incidence rates
case.iNTS$census[year(case.iNTS$date)==2000]<-852054; case.iNTS$census[year(case.iNTS$date)==2001]<-873382
case.iNTS$census[year(case.iNTS$date)==2002]<-894710; case.iNTS$census[year(case.iNTS$date)==2003]<-916039
case.iNTS$census[year(case.iNTS$date)==2004]<-937367; case.iNTS$census[year(case.iNTS$date)==2005]<-958695
case.iNTS$census[year(case.iNTS$date)==2006]<-980023; case.iNTS$census[year(case.iNTS$date)==2007]<-1001352
case.iNTS$census[year(case.iNTS$date)==2008]<-1022680; case.iNTS$census[year(case.iNTS$date)==2009]<-1044008
case.iNTS$census[year(case.iNTS$date)==2010]<-1065337; case.iNTS$census[year(case.iNTS$date)==2011]<-1086665
case.iNTS$census[year(case.iNTS$date)==2012]<-1107993; case.iNTS$census[year(case.iNTS$date)==2013]<-1129322
case.iNTS$census[year(case.iNTS$date)==2014]<-1150650; case.iNTS$census[year(case.iNTS$date)==2015]<-1171978
case.iNTS$incid_sea <-case.iNTS$case_count_sea*100000/case.iNTS$census
case.iNTS$incid_obs <-case.iNTS$case_count_obs*100000/case.iNTS$census
case.typhi$census[year(case.typhi$date)==2000]<-852054; case.typhi$census[year(case.typhi$date)==2001]<-873382
case.typhi$census[year(case.typhi$date)==2002]<-894710; case.typhi$census[year(case.typhi$date)==2003]<-916039
case.typhi$census[year(case.typhi$date)==2004]<-937367; case.typhi$census[year(case.typhi$date)==2005]<-958695
case.typhi$census[year(case.typhi$date)==2006]<-980023; case.typhi$census[year(case.typhi$date)==2007]<-1001352
case.typhi$census[year(case.typhi$date)==2008]<-1022680; case.typhi$census[year(case.typhi$date)==2009]<-1044008
case.typhi$census[year(case.typhi$date)==2010]<-1065337; case.typhi$census[year(case.typhi$date)==2011]<-1086665
case.typhi$census[year(case.typhi$date)==2012]<-1107993; case.typhi$census[year(case.typhi$date)==2013]<-1129322
case.typhi$census[year(case.typhi$date)==2014]<-1150650; case.typhi$census[year(case.typhi$date)==2015]<-1171978
case.typhi$incid_sea <-case.typhi$case_count_sea*100000/case.typhi$census
case.typhi$incid_obs <-case.typhi$case_count_obs*100000/case.typhi$census
Plots of trends in NTS, typhoid, rain and temperature (Supplementary Figure 3)
E1<-ggplot() +
geom_line(data=case.iNTS, aes(x=date, y=case_count_tre, color="Number of NTS cases"), stat="identity", size=1.0) +
geom_line(data=climate.rain, aes(x=date, y=rain_tre/0.04, color="Rainfall"), alpha=0.8, size=0.7) +
scale_y_continuous(sec.axis = sec_axis(~.*0.04, name="Rainfall (mm)")) + theme_bw() +
scale_color_manual(values = c("Number of NTS cases"="orange2", "Rainfall"="blue2")) +
ggtitle("A") + ylab("NTS cases") + xlab("Month'Year") +
theme(axis.title.x = element_text(size=0), axis.title.y=element_text(size=10, color="orange2",face="bold"),
plot.title = element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(size=0), axis.title.y.right=element_text(color="blue2", size=10),
axis.text.y = element_text(face="bold", size=10)) +
scale_x_date(date_breaks = "24 month", date_labels ="%b'%y") +
theme(legend.position="none")
E2<-ggplot() +
geom_line(data=case.iNTS, aes(x=date, y=case_count_tre, color="Number of NTS cases"), stat="identity", size=1.0) +
geom_line(data=climate.temp, aes(x=date, y=temp_tre*4, color="Temperature"), alpha=0.8, size=0.7) +
scale_y_continuous(sec.axis = sec_axis(~./4, name="Temperature (°C)")) + theme_bw() +
scale_color_manual(values=c("Number of NTS cases"="orange2", "Temperature"="green2")) +
ggtitle("B") + ylab("NTS cases") + xlab("Month'Year") +
theme(axis.title.x = element_text(size=10), axis.title.y=element_text(size=10, color="orange2",face="bold"),
plot.title = element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=8), axis.title.y.right=element_text(color="green2", size=10),
axis.text.y=element_text(face="bold", size=10)) +
scale_x_date(date_breaks = "24 month", date_labels ="%b'%y") +
theme(legend.position="none")
E3<-ggplot() +
geom_line(data=case.typhi, aes(x=date, y=case_count_tre, color="Number of typhoid cases"), stat="identity", size=1.0) +
geom_line(data=climate.rain, aes(x=date, y=rain_tre/0.06, color="Rainfall"), alpha=0.8, size=0.7) +
scale_y_continuous(sec.axis=sec_axis(~.*0.06, name="Rainfall (mm)")) + theme_bw() +
scale_color_manual(values=c("Number of typhoid cases"="red2", "Rainfall"="blue2")) +
ggtitle("C") + ylab("Typhoid cases") + xlab("Month'Year") +
theme(axis.title.x = element_text(size=0), axis.title.y=element_text(size=10, color="red2",face="bold"),
plot.title=element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(size=0), axis.title.y.right=element_text(color="blue2", size=10),
axis.text.y = element_text(face="bold", size=10)) +
scale_x_date(date_breaks = "24 month", date_labels ="%b'%y") +
theme(legend.position="none")
E4<-ggplot() +
geom_line(data=case.typhi, aes(x=date, y=case_count_tre, color="Number of typhoid cases"), stat="identity", size=1.0) +
geom_line(data=climate.temp, aes(x=date, y=temp_tre/0.4, color="Temperature"), alpha=0.8, size = 0.7) +
scale_y_continuous(sec.axis=sec_axis(~.*0.4, name="Temperature (°C)")) + theme_bw() +
scale_color_manual(values=c("Number of typhoid cases"="red2", "Temperature"="green2")) +
ggtitle("D") + ylab("Typhoid cases") + xlab("Month'Year") +
theme(axis.title.x=element_text(size=10), axis.title.y=element_text(size=10, color="red2",face="bold"),
plot.title=element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=8), axis.title.y.right = element_text(color="green2", size=10),
axis.text.y = element_text(face="bold", size=10)) +
scale_x_date(date_breaks="24 month", date_labels="%b'%y") +
theme(legend.position="none")
grid.arrange(grobs=list(E1, E3, E2, E4), ncol=2, nrow=2)

Plots of Age-sex distribution of cases (Figure 3)
case$sex[case$sex == ""] <- NA
case$age[case$age == ""] <- NA
case$date <- ymd(case$case_date)
case$year <- year(case$date)
agesex.p1 <- ggplot(subset(case, organism=="iNTS" & !is.na(age) & sex != "Unknown"), aes(x=age, color=sex)) +
geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) +
scale_color_manual(values=c(Male="gray40",Female="red3")) +
theme_bw() +
scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) +
labs(title="A", x="Age (years)", y="Number of NTS cases") +
theme(plot.title=element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) +
theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
labs(color="Missing sex: 2,596 (32.3%)") +
theme(legend.key.height=unit(1,"line")) +
theme(legend.key.width=unit(1,"line"))
agesex.p2<- ggplot(subset(case, organism=="typhi" & !is.na(age) & sex !="Unknown"), aes(x=age, color=sex)) +
geom_freqpoly(position=position_dodge(width=1.5), binwidth=1, size=1) +
scale_color_manual(values=c(Male="gray40",Female="red3")) +
theme_bw() +
scale_x_continuous(limits=c(0, 91), breaks=seq(0, 91, 5)) +
labs(title="B", x="Age (years)", y="Number of typhoid cases") +
theme(plot.title = element_text(hjust=0,face="bold")) +
theme(axis.text.x=element_text(face="bold", size=10, color="black"), axis.text.y=element_text(face="bold", size=10, color="black")) +
theme(legend.justification=c(0.5,0), legend.position=c(0.5, 0.6), legend.text=element_text(size=10), legend.title=element_text(size=10)) +
labs(color="Missing sex: 61 (2.4%)") +
theme(legend.key.height=unit(1,"line")) +
theme(legend.key.width=unit(1,"line"))
grid.arrange(grobs=list(agesex.p1, agesex.p2), ncol=2, nrow=1)

Generating yearly and monthly dynamics
case.iNTS.spi <- subset(case.iNTS, year(case.iNTS$date)<2011, select=c(date,incid_sea))
case.iNTS.spi$month <-month(case.iNTS.spi$date)
case.iNTS.spi$year <- year(case.iNTS.spi$date); case.iNTS.spi$date <- NULL
case.iNTS.spi <- spread(case.iNTS.spi, year, incid_sea); case.iNTS.spi <- as.matrix(case.iNTS.spi)[,-1]
YMD1 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~case.iNTS.spi, type="contour", colorscale='jet', contours=list(showlabels=TRUE)) %>%
colorbar(title="NTS incidence per \n 100,000 population") %>%
layout(title="A", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
case.typhi.spi <- subset(case.typhi, year(case.typhi$date)>2010, select=c(date,incid_sea))
case.typhi.spi$month <- month(case.typhi.spi$date)
case.typhi.spi$year <- year(case.typhi.spi$date); case.typhi.spi$date <- NULL
case.typhi.spi <- spread(case.typhi.spi, year, incid_sea); case.typhi.spi <- as.matrix(case.typhi.spi)[,-1]
YMD2 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~case.typhi.spi, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Typhoid incidence per \n 100,000 population") %>%
layout(title="D", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spin <- subset(climate.rain, year(climate.rain$date)<2011, select=c(date,rainfall_obs))
climate.rain.spin$month <- month(climate.rain.spin$date)
climate.rain.spin$year <- year(climate.rain.spin$date); climate.rain.spin$date <- NULL
climate.rain.spin <- spread(climate.rain.spin, year, rainfall_obs); climate.rain.spin <- as.matrix(climate.rain.spin)[,-1]
YMD3 <- plot_ly(x = c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.rain.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>% colorbar(title="Rainfall (mm)") %>%
layout(title="B", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.rain.spit <- subset(climate.rain, year(climate.rain$date)>2010, select=c(date,rainfall_obs))
climate.rain.spit$month <- month(climate.rain.spit$date)
climate.rain.spit$year <- year(climate.rain.spit$date); climate.rain.spit$date <- NULL
climate.rain.spit <- spread(climate.rain.spit, year, rainfall_obs); climate.rain.spit <- as.matrix(climate.rain.spit)[,-1]
YMD4 <- plot_ly(x = c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.rain.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Rainfall (mm)") %>%
layout(title="E", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spin <- subset(climate.temp, year(climate.temp$date)<2011, select=c(date,temperature_obs))
climate.temp.spin$month <- month(climate.temp.spin$date)
climate.temp.spin$year <- year(climate.temp.spin$date); climate.temp.spin$date <- NULL
climate.temp.spin <- spread(climate.temp.spin, year, temperature_obs); climate.temp.spin <- as.matrix(climate.temp.spin)[,-1]
YMD5 <- plot_ly(x=c(2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.temp.spin, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Temperature (°C)") %>%
layout(title="C", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size=13))
climate.temp.spit <- subset(climate.temp, year(climate.temp$date)>2010, select=c(date,temperature_obs))
climate.temp.spit$month <- month(climate.temp.spit$date)
climate.temp.spit$year <- year(climate.temp.spit$date); climate.temp.spit$date <- NULL
climate.temp.spit <- spread(climate.temp.spit, year, temperature_obs); climate.temp.spit <- as.matrix(climate.temp.spit)[,-1]
YMD6 <- plot_ly(x=c(2011,2012,2013,2014,2015), y=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),
z=~climate.temp.spit, type="contour", colorscale='heatmap', contours=list(showlabels=TRUE)) %>%
colorbar(title="Temperature (°C)") %>%
layout(title="F", xaxis=list(title="Year"), yaxis=list(title="Month"), font=list(size = 13))
Constructing DLNM datasets for NTS and Typhoid
mo.dlnmN <- bind_cols(case.iNTS, climate.rain, climate.temp, id=NULL)
mo.dlnmN$date1 <- mo.dlnmN$date2 <- NULL
mo.dlnmN <- subset(mo.dlnmN, year(date) < 2011)
mo.dlnmN$time <- seq.int(from = 1, to=132, by=1)
mo.dlnmN$year <- year(mo.dlnmN$date)
mo.dlnmN$month <- month(mo.dlnmN$date)
mo.dlnmN$incid_seaX<-round(mo.dlnmN$incid_sea, digits = 0)
mo.dlnmN$incid_obsX<-round(mo.dlnmN$incid_obs, digits = 0)
mo.dlnmT <- bind_cols(case.typhi, climate.rain, climate.temp, id=NULL)
mo.dlnmT$date1 <- mo.dlnmT$date2 <- NULL
mo.dlnmT <- subset(mo.dlnmT, year(date) > 2010)
mo.dlnmT$time <- seq.int(from = 1, to=60, by=1)
mo.dlnmT$year <- year(mo.dlnmT$date)
mo.dlnmT$month <- month(mo.dlnmT$date)
mo.dlnmT$incid_seaX<-round(mo.dlnmT$incid_sea, digits = 0)
Switching from AIC to QAIC for NTS and Typhoid
nts.quasipoisson <- function(...) {
res <- quasipoisson(...)
res$aic <- poisson(...)$aic
res
}
typ.quasipoisson <- function(...) {
res <- quasipoisson(...)
res$aic <- poisson(...)$aic
res
}
Calculate QAIC values to select optimal models (Supplementary Table)
QAICtable <- data.frame(model.no=rep(NA,27), lag.df=rep(NA,27), fx1.df=rep(NA,27), fx2.df=rep(NA,27),
QAIC.ntsR=rep(NA,27), QAIC.ntsT=rep(NA,27), QAIC.nts=rep(NA,27), QAIC.typR=rep(NA,27),
QAIC.typT=rep(NA,27), QAIC.typ=rep(NA,27))
l=1
for(k in 3:5){
for(j in 3:5){
for(i in 3:5){
nts.lagknots <- logknots(8, fun="ns", df=i)
nts.varknots.r=equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=j)
nts.varknots.t=equalknots(mo.dlnmN$temperature_obs, fun="ns", df=k)
nts.mo.cb.rain <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots))
nts.mo.cb.temp <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots))
nts.modelR <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
nts.modelT <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.temp + year, family=nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
nts.model <- glm(mo.dlnmN$incid_seaX ~ nts.mo.cb.rain + nts.mo.cb.temp + year, family = nts.quasipoisson(), na.action=na.delete, mo.dlnmN)
typ.lagknots <- logknots(8, fun="ns", df=i)
typ.varknots.r=equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=j)
typ.varknots.t=equalknots(mo.dlnmT$temperature_obs, fun="ns", df=k)
typ.mo.cb.rain <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r), arglag=list(knots=typ.lagknots))
typ.mo.cb.temp <- crossbasis(mo.dlnmT$temperature_obs, lag =8, argvar = list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.model <- glm(mo.dlnmT$incid_seaX ~ typ.mo.cb.rain + typ.mo.cb.temp + year, family = typ.quasipoisson(), na.action=na.delete, mo.dlnmT)
QAICtable[l,] <- c(l,i,j,k, QAIC(nts.modelR, chat=summary(nts.modelR)$dispersion), QAIC(nts.modelT, chat=summary(nts.modelT)$dispersion),
QAIC(nts.model, chat=summary(nts.model)$dispersion), QAIC(typ.modelR, chat=summary(typ.modelR)$dispersion), QAIC(typ.modelT, chat=summary(typ.modelT)$dispersion),
QAIC(typ.model, chat=1))
l=l+1
}
}
}
kable(QAICtable)
| 1 |
3 |
3 |
3 |
570.8194 |
565.8631 |
562.3950 |
218.5004 |
220.2006 |
212.6658 |
| 2 |
4 |
3 |
3 |
576.3831 |
571.4838 |
573.3630 |
223.1714 |
216.1853 |
220.6390 |
| 3 |
5 |
3 |
3 |
581.4852 |
574.9872 |
583.4617 |
228.2399 |
221.1698 |
230.8468 |
| 4 |
3 |
4 |
3 |
571.1280 |
565.8631 |
565.3731 |
211.3436 |
220.2006 |
216.2219 |
| 5 |
4 |
4 |
3 |
578.3398 |
571.4838 |
578.1025 |
218.1843 |
216.1853 |
228.0142 |
| 6 |
5 |
4 |
3 |
585.4681 |
574.9872 |
590.3358 |
222.7394 |
221.1698 |
237.6362 |
| 7 |
3 |
5 |
3 |
572.7386 |
565.8631 |
566.1964 |
216.3176 |
220.2006 |
221.5570 |
| 8 |
4 |
5 |
3 |
581.6967 |
571.4838 |
581.3661 |
222.4327 |
216.1853 |
234.0871 |
| 9 |
5 |
5 |
3 |
590.5606 |
574.9872 |
594.8792 |
227.8628 |
221.1698 |
246.2023 |
| 10 |
3 |
3 |
4 |
570.8194 |
571.7243 |
566.2250 |
218.5004 |
220.7893 |
215.5219 |
| 11 |
4 |
3 |
4 |
576.3831 |
577.5401 |
578.9833 |
223.1714 |
220.2852 |
225.9024 |
| 12 |
5 |
3 |
4 |
581.4852 |
583.4447 |
589.9669 |
228.2399 |
226.2618 |
236.0380 |
| 13 |
3 |
4 |
4 |
571.1280 |
571.7243 |
568.6740 |
211.3436 |
220.7893 |
220.0702 |
| 14 |
4 |
4 |
4 |
578.3398 |
577.5401 |
583.3435 |
218.1843 |
220.2852 |
233.1946 |
| 15 |
5 |
4 |
4 |
585.4681 |
583.4447 |
597.5272 |
222.7394 |
226.2618 |
245.6773 |
| 16 |
3 |
5 |
4 |
572.7386 |
571.7243 |
571.2271 |
216.3176 |
220.7893 |
224.5520 |
| 17 |
4 |
5 |
4 |
581.6967 |
577.5401 |
588.2408 |
222.4327 |
220.2852 |
239.9495 |
| 18 |
5 |
5 |
4 |
590.5606 |
583.4447 |
603.1850 |
227.8628 |
226.2618 |
252.6551 |
| 19 |
3 |
3 |
5 |
570.8194 |
574.0369 |
566.7772 |
218.5004 |
220.7423 |
220.4298 |
| 20 |
4 |
3 |
5 |
576.3831 |
580.7484 |
580.7418 |
223.1714 |
220.9285 |
232.7209 |
| 21 |
5 |
3 |
5 |
581.4852 |
587.7540 |
593.8723 |
228.2399 |
225.5709 |
245.1746 |
| 22 |
3 |
4 |
5 |
571.1280 |
574.0369 |
569.9469 |
211.3436 |
220.7423 |
224.8192 |
| 23 |
4 |
4 |
5 |
578.3398 |
580.7484 |
586.2703 |
218.1843 |
220.9285 |
239.8864 |
| 24 |
5 |
4 |
5 |
585.4681 |
587.7540 |
601.6447 |
222.7394 |
225.5709 |
254.6679 |
| 25 |
3 |
5 |
5 |
572.7386 |
574.0369 |
575.0626 |
216.3176 |
220.7423 |
229.8974 |
| 26 |
4 |
5 |
5 |
581.6967 |
580.7484 |
593.7784 |
222.4327 |
220.9285 |
246.7405 |
| 27 |
5 |
5 |
5 |
590.5606 |
587.7540 |
610.6531 |
227.8628 |
225.5709 |
261.6885 |
Fit NTS model
nts.model <- glm(mo.dlnmN$incid_seaX ~ nts.cb.rain + nts.cb.temp + year, family=quasipoisson(), na.action=na.delete, mo.dlnmN)
Investigate NTS deviance residuals
par(mfrow=c(1,2))
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Original model with Rain & Temp", xlim=c(0,8))
nts.model <- update(nts.model,.~.+Lag(residuals(nts.model,type="deviance"), 1)) #add residuals at lag 1
pacf(residuals(nts.model,type="deviance"), na.action=na.omit, main="Adjusted model with Rain & Temp", xlim=c(0,8))

Run NTS predictions
nts.pred.rain <- crosspred(nts.cb.rain, nts.model, cen=0, by=0.2)
nts.pred.temp <- crosspred(nts.cb.temp, nts.model, cen=23, by=0.2)
Plots of NTS predictions due to rainfall (Figure 4)
#plot(nts.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1.1, cex.axis=1.1, col="gray80",main="A")
plot(nts.pred.rain, "contour", key.title=title("NTS"), plot.title=title("", xlab ="Rainfall (mm)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="A"))

plot(nts.pred.rain, "slices", xlab="Lag (month) [given 9 mm] ", var=c(9), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="B")

plot(nts.pred.rain, "slices", xlab="Lag (month) [given 13 mm] ", var=c(13), col="orange2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95, ci='b', lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="C")

Plots of NTS predictions due to temperature (Figure 6)
#plot(nts.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of NTS", theta=150, phi=5, lphi=100, cex.lab=1, cex.axis=1, col="gray80",main="A")
plot(nts.pred.temp, "contour", key.title=title("NTS"), plot.title=title("", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="A"))

plot(nts.pred.temp, xlab="Lag (month) [given 19 °C]", "slices",var=c(19), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="B")

plot(nts.pred.temp, xlab="Lag (month) [given 29 °C]", "slices",var=c(29), col="orange2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95, ci='b',lwd=4.5, ylab="RR of NTS", cex.lab=1.5, cex.axis=1.5,main="C")

Fit Typhoid models
typ.modelR <- glm(mo.dlnmT$incid_seaX ~ typ.cb.rain1 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
typ.modelT <- glm(mo.dlnmT$incid_seaX ~ typ.cb.temp + typ.cb.rain2 + year, family=quasipoisson(), na.action=na.delete, mo.dlnmT)
Investigate deviance residuals for fitted tyhpoid models
par(mfrow=c(1,3))
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="Original model with Rain",xlim=c(0,8))
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Original model with Rain & Temp",xlim=c(0,8))
typ.modelT <- update(typ.modelT,.~.+Lag(residuals(typ.modelT,type="deviance"),1)) #add residuals at lag 1
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="Adjusted model with Rain & Temp",xlim=c(0,8))

Run typhoid predictions
typ.pred.rain <- crosspred(typ.cb.rain1, typ.modelR, cen=0, by=0.2)
typ.pred.temp <- crosspred(typ.cb.temp, typ.modelT, cen=23, by=0.2)
Plots of TYP predictions due to rainfall (Figure 5)
#plot(typ.pred.rain, xlab="Rainfall (mm)", ylab="Lag (month)", zlab="RR of typhoid", theta=210, phi=10, lphi=150, col="gray80", main="D", cex.lab=1.1, cex.axis=1.1)
plot(typ.pred.rain, "contour", key.title=title("TYP"), plot.title=title("", xlab ="Rainfall (mm)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="D"))

plot(typ.pred.rain, xlab="Lag (month) [given 9 mm]", "slices", var=c(9), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="E")

plot(typ.pred.rain, xlab="Lag (month) [given 16 mm]", "slices", var=c(16), col="red2", ci.arg=list(col=topo.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="F")

Plots of TYP predictions due to temperature (Figure 7)
#plot(typ.pred.temp, xlab="Temperature (°C)", ylab="Lag (month)", zlab="RR of typhoid", theta=150, phi=5, lphi=100, col="gray80", main="D", cex.lab=1.1, cex.axis=1.1)
plot(typ.pred.temp, "contour", key.title=title("TYP"), plot.title=title("", xlab ="Temperature (°C)", ylab = "Lag (month)", cex.lab=1.7, cex.axis=1.7,main="D"))

plot(typ.pred.temp, xlab="Lag (month) [given 19 °C]", "slices", var=c(19), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="E")

plot(typ.pred.temp, xlab="Lag (month) [given 25 °C]", "slices", var=c(25), col="red2", ci.arg=list(col=terrain.colors(70, alpha = 1)), ci.level=0.95,ci='b', lwd=4.5, ylab="RR of typhoid", cex.lab=1.5, cex.axis=1.5,main="F")

Plots of validation checks (Supplementary Figure S4)
par(mfrow=c(3,4))
plot(mo.dlnmN$date, residuals(nts.model,type="deviance"), pch=19, cex=0.8, col=grey(0.4),main="A", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2)
abline(h=0,lty=2,lwd=3)
pacf(residuals(nts.model,type="deviance"),na.action=na.omit,main="B",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
acf(residuals(nts.model,type="deviance"),na.action=na.omit,main="C",xlim=c(0,8),xlab="", cex.lab=1.2,cex.axis=1.2)
plot(matrix(mo.dlnmN$incid_seaX),matrix(predict(nts.model, type="response")), col=c("orange2","gray2"),main="D",xlab="Predicted NTS incidence", ylab="Observed NTS incidence",pch=19,cex=0.8,cex.lab=1.2,cex.axis=1.2)
legend("topleft", legend=c("observed", "predicted"), col=c("grey30", "orange2"), cex=0.8, pch=19)
plot(mo.dlnmT$date, residuals(typ.modelR,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="E", ylab="Residuals",xlab="",cex.lab=1.2,cex.axis=1.2)
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="F",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelR,type="deviance"),na.action=na.omit,main="G",xlim=c(0,8),xlab="",cex.lab=1.2, cex.axis=1.2)
plot(matrix(mo.dlnmT$incid_seaX),matrix(predict(typ.modelR, type="response")), col=c("red2","gray2"),main="H",xlab="Predicted typhoid incidence", ylab="Observed typhoid incidence",pch=19,cex=0.8,cex.lab=1.2,cex.axis=1.2)
legend("topleft", legend=c("observed", "predicted"), col=c("grey30", "red2"), cex=0.8, pch=19)
plot(mo.dlnmT$date, residuals(typ.modelT,type="deviance"), pch=19, cex=0.8, col=grey(0.6),main="I", ylab="Residuals", xlab="Year",cex.lab=1.2,cex.axis=1.2)
abline(h=0,lty=2,lwd=3)
pacf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="J",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
acf(residuals(typ.modelT,type="deviance"),na.action=na.omit,main="K",xlim=c(0,8),xlab="Lag (month)",cex.lab=1.2, cex.axis=1.2)
plot(matrix(mo.dlnmT$incid_seaX),matrix(predict(typ.modelT, type="response")), col=c("red2","gray2"),main="L",xlab="Predicted typhoid incidence", ylab="Observed typhoid incidence",pch=19,cex=0.8,cex.lab=1.2,cex.axis=1.2)
legend("topleft", legend=c("observed", "predicted"), col=c("grey30", "red2"), cex=0.8, pch=19)

Plots of predictions due to rainfall/temperature (Table 1)
RRTable <- data.frame(rbind(nts.pred.rain$matRRfit["9",], nts.pred.rain$matRRlow["9",], nts.pred.rain$matRRhigh["9",],
nts.pred.rain$matRRfit["13",],nts.pred.rain$matRRlow["13",],nts.pred.rain$matRRhigh["13",],
nts.pred.temp$matRRfit["19",],nts.pred.temp$matRRlow["19",],nts.pred.temp$matRRhigh["19",],
nts.pred.temp$matRRfit["29",],nts.pred.temp$matRRlow["29",],nts.pred.temp$matRRhigh["29",],
typ.pred.rain$matRRfit["9",],typ.pred.rain$matRRlow["9",],typ.pred.rain$matRRhigh["9",],
typ.pred.rain$matRRfit["16",],typ.pred.rain$matRRlow["16",],typ.pred.rain$matRRhigh["16",],
typ.pred.temp$matRRfit["19",],typ.pred.temp$matRRlow["19",],typ.pred.temp$matRRhigh["19",],
typ.pred.temp$matRRfit["25",],typ.pred.temp$matRRlow["25",],typ.pred.temp$matRRhigh["25",]
))
kable(RRTable)
| 1.2248467 |
1.2514626 |
1.2554780 |
1.2238951 |
1.1635420 |
1.0841819 |
0.9951394 |
0.9042881 |
0.8176177 |
| 1.0254867 |
1.1162267 |
1.1453653 |
1.1143857 |
1.0596817 |
0.9906237 |
0.9027750 |
0.7975713 |
0.6897005 |
| 1.4629635 |
1.4030831 |
1.3761766 |
1.3441658 |
1.2775817 |
1.1865762 |
1.0969538 |
1.0252839 |
0.9692593 |
| 1.1817325 |
1.0655402 |
0.9723354 |
0.9040825 |
0.8545278 |
0.8183618 |
0.7914806 |
0.7705230 |
0.7525862 |
| 0.9099591 |
0.8927064 |
0.8190710 |
0.7501790 |
0.7124350 |
0.6968964 |
0.6814001 |
0.6417570 |
0.5822991 |
| 1.5346752 |
1.2718359 |
1.1542787 |
1.0895601 |
1.0249607 |
0.9609980 |
0.9193445 |
0.9251255 |
0.9726719 |
| 1.0822001 |
1.1309828 |
1.1629884 |
1.1659357 |
1.1432279 |
1.1012306 |
1.0467408 |
0.9861515 |
0.9249535 |
| 0.9602616 |
1.0318836 |
1.0674043 |
1.0690932 |
1.0514755 |
1.0192753 |
0.9708904 |
0.9049599 |
0.8287981 |
| 1.2196229 |
1.2395991 |
1.2671318 |
1.2715505 |
1.2429867 |
1.1897756 |
1.1285169 |
1.0746274 |
1.0322647 |
| 0.7825270 |
0.7912493 |
0.7912237 |
0.7775272 |
0.7525056 |
0.7194590 |
0.6815991 |
0.6418050 |
0.6024941 |
| 0.6253076 |
0.6887091 |
0.6928915 |
0.6692501 |
0.6487554 |
0.6292260 |
0.5946514 |
0.5350856 |
0.4650197 |
| 0.9792756 |
0.9090565 |
0.9035107 |
0.9033223 |
0.8728476 |
0.8226317 |
0.7812600 |
0.7698090 |
0.7806102 |
| 0.9870340 |
1.2979004 |
1.6063110 |
1.8078055 |
1.8722857 |
1.8143089 |
1.6725894 |
1.4915122 |
1.3081090 |
| 0.6284634 |
0.9438326 |
1.2552773 |
1.4434731 |
1.5069006 |
1.4634299 |
1.3249680 |
1.1234975 |
0.9124989 |
| 1.5501876 |
1.7847926 |
2.0555100 |
2.2640954 |
2.3262675 |
2.2493164 |
2.1114137 |
1.9800744 |
1.8752341 |
| 0.2199818 |
0.2071173 |
0.2101761 |
0.2398628 |
0.3033687 |
0.4165633 |
0.6083705 |
0.9257774 |
1.4380388 |
| 0.0679390 |
0.0908122 |
0.1004582 |
0.1120403 |
0.1444326 |
0.2128310 |
0.3269832 |
0.4616799 |
0.5759679 |
| 0.7122860 |
0.4723769 |
0.4397253 |
0.5135131 |
0.6372006 |
0.8153182 |
1.1319073 |
1.8564025 |
3.5904008 |
| 1.4960476 |
1.5619360 |
1.5874626 |
1.5468052 |
1.4526163 |
1.3245023 |
1.1812571 |
1.0380766 |
0.9055470 |
| 1.0424914 |
1.1457962 |
1.0968536 |
1.0279449 |
0.9761052 |
0.9252186 |
0.8373032 |
0.6896236 |
0.5229565 |
| 2.1469324 |
2.1292129 |
2.2975149 |
2.3275628 |
2.1617487 |
1.8960993 |
1.6665030 |
1.5625959 |
1.5680377 |
| 0.6881750 |
0.9207072 |
1.1640344 |
1.3467161 |
1.4416997 |
1.4504554 |
1.3928588 |
1.2966573 |
1.1885051 |
| 0.4844865 |
0.7123388 |
0.8891649 |
0.9975957 |
1.0683040 |
1.1081654 |
1.1070709 |
1.0400232 |
0.9069586 |
| 0.9774984 |
1.1900260 |
1.5238750 |
1.8180153 |
1.9456054 |
1.8984720 |
1.7524222 |
1.6166181 |
1.5574519 |
Vary degrees of freedom for sensitivity analyses of NTS/typhoid predictions
nts.cb.rain <- list()
nts.cb.temp <- list()
typ.cb.rain <- list()
typ.cb.temp <- list()
sort(mo.dlnmN$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmN$temperature_obs, decreasing=FALSE)
sort(mo.dlnmT$rainfall_obs, decreasing=FALSE)
sort(mo.dlnmT$temperature_obs, decreasing=FALSE)
l=1
for(k in 3:5){
for(j in 3:5){
for(i in 3:5){
nts.lagknots <- logknots(8, fun="ns", df=i)
nts.varknots.r=equalknots(mo.dlnmN$rainfall_obs, fun="ns", df=j)
nts.varknots.t=equalknots(mo.dlnmN$temperature_obs, fun="ns", df=k)
nts.cb.rain[[l]] <- crossbasis(mo.dlnmN$rainfall_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.r), arglag=list(knots=nts.lagknots))
nts.cb.temp[[l]] <- crossbasis(mo.dlnmN$temperature_obs, lag=8, argvar=list(fun="ns", knots=nts.varknots.t), arglag=list(knots=nts.lagknots))
typ.lagknots <- logknots(8, fun="ns", df=i)
typ.varknots.r=equalknots(mo.dlnmT$rainfall_obs, fun="ns", df=j)
typ.varknots.t=equalknots(mo.dlnmT$temperature_obs, fun="ns", df=k)
typ.cb.rain[[l]] <- crossbasis(mo.dlnmT$rainfall_obs, lag=8, argvar=list(fun="ns", knots=typ.varknots.r), arglag=list(knots=typ.lagknots))
typ.cb.temp[[l]] <- crossbasis(mo.dlnmT$temperature_obs, lag =8, argvar = list(fun="ns", knots=typ.varknots.t), arglag=list(knots=typ.lagknots))
l=l+1
}
}
}
Formulate alternative sensitivity models and plot predictions of NTS (Supplementay Figure S5)
nts.cb.rain.s1 <- nts.cb.rain[[2]]; nts.cb.temp.s1 <- nts.cb.temp[[2]]
nts.cb.rain.s2 <- nts.cb.rain[[3]]; nts.cb.temp.s2 <- nts.cb.temp[[3]]
nts.cb.rain.s3 <- nts.cb.rain[[4]]; nts.cb.temp.s3 <- nts.cb.temp[[4]]
nts.cb.rain.s4 <- nts.cb.rain[[7]]; nts.cb.temp.s4 <- nts.cb.temp[[7]]
nts.cb.rain.s5 <- nts.cb.rain[[10]]; nts.cb.temp.s5 <- nts.cb.temp[[10]]
nts.cb.rain.s6 <- nts.cb.rain[[19]]; nts.cb.temp.s6 <- nts.cb.temp[[19]]
nts.model.s1 <- glm(incid_seaX~nts.cb.rain.s1 + nts.cb.temp.s1 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s2 <- glm(incid_seaX~nts.cb.rain.s2 + nts.cb.temp.s2 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s3 <- glm(incid_seaX~nts.cb.rain.s3 + nts.cb.temp.s3 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s4 <- glm(incid_seaX~nts.cb.rain.s4 + nts.cb.temp.s4 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s5 <- glm(incid_seaX~nts.cb.rain.s5 + nts.cb.temp.s5 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s6 <- glm(incid_seaX~nts.cb.rain.s6 + nts.cb.temp.s6 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmN)
nts.model.s1 <- update(nts.model.s1,.~.+Lag(residuals(nts.model.s1,type="deviance"),1)) #add residuals at lag 1
nts.model.s2 <- update(nts.model.s2,.~.+Lag(residuals(nts.model.s2,type="deviance"),1)) #to correct for partial
nts.model.s3 <- update(nts.model.s3,.~.+Lag(residuals(nts.model.s3,type="deviance"),1)) #autocorrelation
nts.model.s4 <- update(nts.model.s4,.~.+Lag(residuals(nts.model.s4,type="deviance"),1))
nts.model.s5 <- update(nts.model.s5,.~.+Lag(residuals(nts.model.s5,type="deviance"),1))
nts.model.s6 <- update(nts.model.s6,.~.+Lag(residuals(nts.model.s6,type="deviance"),1))
nts.pred.rain.s1 <- crosspred(nts.cb.rain.s1, nts.model.s1, cen=0, by=0.2)
nts.pred.rain.s2 <- crosspred(nts.cb.rain.s2, nts.model.s2, cen=0, by=0.2)
nts.pred.rain.s3 <- crosspred(nts.cb.rain.s3, nts.model.s3, cen=0, by=0.2)
nts.pred.rain.s4 <- crosspred(nts.cb.rain.s4, nts.model.s4, cen=0, by=0.2)
nts.pred.rain.s5 <- crosspred(nts.cb.rain.s5, nts.model.s5, cen=0, by=0.2)
nts.pred.rain.s6 <- crosspred(nts.cb.rain.s6, nts.model.s6, cen=0, by=0.2)
nts.pred.temp.s1 <- crosspred(nts.cb.temp.s1, nts.model.s1, cen=23, by=0.2)
nts.pred.temp.s2 <- crosspred(nts.cb.temp.s2, nts.model.s2, cen=23, by=0.2)
nts.pred.temp.s3 <- crosspred(nts.cb.temp.s3, nts.model.s3, cen=23, by=0.2)
nts.pred.temp.s4 <- crosspred(nts.cb.temp.s4, nts.model.s4, cen=23, by=0.2)
nts.pred.temp.s5 <- crosspred(nts.cb.temp.s5, nts.model.s5, cen=23, by=0.2)
nts.pred.temp.s6 <- crosspred(nts.cb.temp.s6, nts.model.s6, cen=23, by=0.2)
par(mfrow=c(2,6))
plot(nts.pred.rain.s1,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="A")
plot(nts.pred.rain.s2,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="C")
plot(nts.pred.rain.s3,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="E")
plot(nts.pred.rain.s4,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="G")
plot(nts.pred.rain.s5,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="I")
plot(nts.pred.rain.s6,"slices",var=c(9),col="orange2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="K")
plot(nts.pred.temp.s1,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="B")
plot(nts.pred.temp.s2,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="D")
plot(nts.pred.temp.s3,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="F")
plot(nts.pred.temp.s4,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="H")
plot(nts.pred.temp.s5,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="J")
plot(nts.pred.temp.s6,"slices",var=c(19),col="orange2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of NTS",cex.lab=1.1,cex.axis=1.1,main="L")

Formulate alternative sensitivity models and plot predictions of Typhoid (Supplementay Figure S6)
typ.cb.rain.s1 <- typ.cb.rain[[2]]; typ.cb.temp.s1 <- typ.cb.temp[[2]]
typ.cb.rain.s2 <- typ.cb.rain[[3]]; typ.cb.temp.s2 <- typ.cb.temp[[3]]
typ.cb.rain.s3 <- typ.cb.rain[[4]]; typ.cb.temp.s3 <- typ.cb.temp[[4]]
typ.cb.rain.s4 <- typ.cb.rain[[7]]; typ.cb.temp.s4 <- typ.cb.temp[[7]]
typ.cb.rain.s5 <- typ.cb.rain[[10]]; typ.cb.temp.s5 <- typ.cb.temp[[10]]
typ.cb.rain.s6 <- typ.cb.rain[[19]]; typ.cb.temp.s6 <- typ.cb.temp[[19]]
typ.modelR.s1 <- glm(incid_seaX~typ.cb.rain.s1 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s2 <- glm(incid_seaX~typ.cb.rain.s2 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s3 <- glm(incid_seaX~typ.cb.rain.s3 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s4 <- glm(incid_seaX~typ.cb.rain.s4 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s5 <- glm(incid_seaX~typ.cb.rain.s5 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s6 <- glm(incid_seaX~typ.cb.rain.s6 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s1 <- glm(incid_seaX~typ.cb.rain.s1 + typ.cb.temp.s1 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s2 <- glm(incid_seaX~typ.cb.rain.s2 + typ.cb.temp.s2 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s3 <- glm(incid_seaX~typ.cb.rain.s3 + typ.cb.temp.s3 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s4 <- glm(incid_seaX~typ.cb.rain.s4 + typ.cb.temp.s4 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s5 <- glm(incid_seaX~typ.cb.rain.s5 + typ.cb.temp.s5 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelT.s6 <- glm(incid_seaX~typ.cb.rain.s6 + typ.cb.temp.s6 + year, family=quasipoisson(), na.action=na.exclude, mo.dlnmT)
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),1))
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),2))
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),3))
typ.modelR.s1 <- update(typ.modelR.s1,.~.+Lag(residuals(typ.modelR.s1,type="deviance"),4))
typ.modelR.s2 <- update(typ.modelR.s2,.~.+Lag(residuals(typ.modelR.s2,type="deviance"),1))
typ.modelR.s2 <- update(typ.modelR.s2,.~.+Lag(residuals(typ.modelR.s2,type="deviance"),2))
typ.modelR.s2 <- update(typ.modelR.s2,.~.+Lag(residuals(typ.modelR.s2,type="deviance"),3))
typ.modelR.s5 <- update(typ.modelR.s5,.~.+Lag(residuals(typ.modelR.s5,type="deviance"),1))
typ.modelR.s5 <- update(typ.modelR.s5,.~.+Lag(residuals(typ.modelR.s5,type="deviance"),2))
typ.modelR.s5 <- update(typ.modelR.s5,.~.+Lag(residuals(typ.modelR.s5,type="deviance"),3))
typ.modelR.s6 <- update(typ.modelR.s6,.~.+Lag(residuals(typ.modelR.s6,type="deviance"),1))
typ.modelR.s6 <- update(typ.modelR.s6,.~.+Lag(residuals(typ.modelR.s6,type="deviance"),2))
typ.modelR.s6 <- update(typ.modelR.s6,.~.+Lag(residuals(typ.modelR.s6,type="deviance"),3))
typ.modelT.s1 <- update(typ.modelT.s1,.~.+Lag(residuals(typ.modelT.s1,type="deviance"),1))
typ.modelT.s1 <- update(typ.modelT.s1,.~.+Lag(residuals(typ.modelT.s1,type="deviance"),3))
typ.modelT.s2 <- update(typ.modelT.s2,.~.+Lag(residuals(typ.modelT.s2,type="deviance"),1))
typ.modelT.s3 <- update(typ.modelT.s3,.~.+Lag(residuals(typ.modelT.s3,type="deviance"),1))
typ.modelT.s4 <- update(typ.modelT.s4,.~.+Lag(residuals(typ.modelT.s4,type="deviance"),1))
typ.modelT.s4 <- update(typ.modelT.s4,.~.+Lag(residuals(typ.modelT.s4,type="deviance"),2))
typ.modelT.s5 <- update(typ.modelT.s5,.~.+Lag(residuals(typ.modelT.s5,type="deviance"),1))
typ.modelT.s5 <- update(typ.modelT.s5,.~.+Lag(residuals(typ.modelT.s5,type="deviance"),4))
typ.modelT.s5 <- update(typ.modelT.s5,.~.+Lag(residuals(typ.modelT.s5,type="deviance"),3))
typ.modelT.s6 <- update(typ.modelT.s6,.~.+Lag(residuals(typ.modelT.s6,type="deviance"),1))
typ.pred.rain.s1 <- crosspred(typ.cb.rain.s1, typ.modelR.s1, cen=0, by=0.2)
typ.pred.rain.s2 <- crosspred(typ.cb.rain.s2, typ.modelR.s2, cen=0, by=0.2)
typ.pred.rain.s3 <- crosspred(typ.cb.rain.s3, typ.modelR.s3, cen=0, by=0.2)
typ.pred.rain.s4 <- crosspred(typ.cb.rain.s4, typ.modelR.s4, cen=0, by=0.2)
typ.pred.rain.s5 <- crosspred(typ.cb.rain.s5, typ.modelR.s5, cen=0, by=0.2)
typ.pred.rain.s6 <- crosspred(typ.cb.rain.s6, typ.modelR.s6, cen=0, by=0.2)
typ.pred.temp.s1 <- crosspred(typ.cb.temp.s1, typ.modelT.s1, cen=23, by=0.2)
typ.pred.temp.s2 <- crosspred(typ.cb.temp.s2, typ.modelT.s2, cen=23, by=0.2)
typ.pred.temp.s3 <- crosspred(typ.cb.temp.s3, typ.modelT.s3, cen=23, by=0.2)
typ.pred.temp.s4 <- crosspred(typ.cb.temp.s4, typ.modelT.s4, cen=23, by=0.2)
typ.pred.temp.s5 <- crosspred(typ.cb.temp.s5, typ.modelT.s5, cen=23, by=0.2)
typ.pred.temp.s6 <- crosspred(typ.cb.temp.s6, typ.modelT.s6, cen=23, by=0.2)
par(mfrow=c(2,6))
plot(typ.pred.rain.s1,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="A")
plot(typ.pred.rain.s2,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="C")
plot(typ.pred.rain.s3,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="E")
plot(typ.pred.rain.s4,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="G")
plot(typ.pred.rain.s5,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="I")
plot(typ.pred.rain.s6,"slices",var=c(9),col="red2",ci.arg=list(col=topo.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 9 mm)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="K")
plot(typ.pred.temp.s1,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR oftyphoid",cex.lab=1.1,cex.axis=1.1,main="B")
plot(typ.pred.temp.s2,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="D")
plot(typ.pred.temp.s3,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="F")
plot(typ.pred.temp.s4,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="H")
plot(typ.pred.temp.s5,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="J")
plot(typ.pred.temp.s6,"slices",var=c(19),col="red2",ci.arg=list(col=terrain.colors(70,alpha=1)),ci.level=0.95,ci='b',lwd=4.5,xlab="Month-lag (given 19 °C)",ylab="RR of typhoid",cex.lab=1.1,cex.axis=1.1,main="L")

Descriptive stats of the study population
mean(mo.dlnmN$incid_seaX)
## [1] 6.795455
ci(mo.dlnmN$incid_seaX)
## Estimate CI lower CI upper Std. Error
## 6.7954545 6.1608513 7.4300578 0.3207921
mean(mo.dlnmT$incid_seaX)
## [1] 3.633333
ci(mo.dlnmT$incid_seaX)
## Estimate CI lower CI upper Std. Error
## 3.6333333 2.9011782 4.3654884 0.3658954